The evolutionary history of Brachyury genes in Hydrozoa involves duplications, divergence, and neofunctionalization

Brachyury, a member of T-box gene family, is widely known for its major role in mesoderm specification in bilaterians. It is also present in non-bilaterian metazoans, such as cnidarians, where it acts as a component of an axial patterning system. In this study, we present a phylogenetic analysis of Brachyury genes within phylum Cnidaria, investigate differential expression and address a functional framework of Brachyury paralogs in hydrozoan Dynamena pumila. Our analysis indicates two duplication events of Brachyury within the cnidarian lineage. The first duplication likely appeared in the medusozoan ancestor, resulting in two copies in medusozoans, while the second duplication arose in the hydrozoan ancestor, resulting in three copies in hydrozoans. Brachyury1 and 2 display a conservative expression pattern marking the oral pole of the body axis in D. pumila. On the contrary, Brachyury3 expression was detected in scattered presumably nerve cells of the D. pumila larva. Pharmacological modulations indicated that Brachyury3 is not under regulation of cWnt signaling in contrast to the other two Brachyury genes. Divergence in expression patterns and regulation suggest neofunctionalization of Brachyury3 in hydrozoans.


Results
Diversity and phylogeny of cnidarian brachyury genes. To address the evolution of the Brachyury gene family within cnidarians, first we conducted TBLASTX search of the previously published transcriptome of D. pumila 35 with the published C. hemisphaerica Brachyury gene sequences as an initial query. We recovered three sequences of Brachyury-like genes from the D. pumila transcriptome and used them as queries for TBLASTX searches against ten more medusozoan transcriptomes (see "Methods"). Together with four already known anthozoan sequences, a total of 33 Brachyury sequences were identified from 16 Cnidaria species.
A maximum likelihood tree was generated using translated amino acid sequences with the best-fit JTT++ R5 model (Fig. 1). For this analysis, a total of 41 Brachyury sequences were used representing all major metazoan groups except Porifera. Since the T-box transcription factor family includes classes of Tbx genes besides Brachyury 9 , sequences of metazoan Tbx genes were used as an out-group to root the tree. To test the robustness of the tree topology, we also used the conservative T-box domain for the alignment. Additional maximum likelihood tree with the same overall topology was generated using only T-boxes of analyzed sequences (Supplementary Information Fig. S1).
For a more thorough analysis of Brachyury genes in Hydra we searched for Brachyury transcripts in gene models of Hydra 2.0 (Hydra magnipapillata) 37 and HydraAEP (Hydra vulgaris) 38 genome assemblies with modified PIA3 pipeline, which allowed us to retrieve T-box protein class information automatically. This analysis confirmed that Hydra genomes contain only two Brachyury genes and several T-box genes (see an example of tree output in Supplementary Information Fig. S2). The manual search also retrieved only two Brachyury genes in Hydra. All these results increase the likelihood that Hydra indeed lost Brachyury3.
Comparison of sequence conservation of hydrozoan Brachyury proteins. Multiple sequence alignment of the deduced full-length amino acid sequences of D. pumila Brachyury proteins revealed that their T-boxes show about 70-77% identity. By contrast, full-length sequences have an overall lower amino acid identity, thus, the remaining regions are less conserved (Fig. 2a, b).
Further, we performed a multiple sequence alignment of deduced Brachyury proteins, functional domain prediction by hmmscan, and a search of the R1 repressor domain with a sequence of R1 from H. vulgaris Brachy-ury1 as query 33 . Hydrozoan Brachyury1 proteins share higher identities with homologous genes in different species than Brachyury2 and Brachyury3 (Fig. 2c, d). The inter-species Brachyury1 identities were also higher than identities between Brachyury paralogs within the same species (Fig. 2b, d). Only Brachyury1 proteins have the repression R1 fragment in the C-terminal region (Fig. 2c, Supplementary Information Fig. S3). Scyphozoan and cubozoan Brachyury2/3 proteins ( Supplementary Information Fig. S3) and hydrozoan Brachyury2 and Brachyury3 proteins (Fig. 2c) have lost it. Inter-species sequence comparison further revealed that of all analyzed hydrozoan Brachyury proteins, Brachyury2 proteins have the shortest sequences positioned N-terminally of the T-box (e.g., 20 amino acid for DpBra2), whereas Brachyury3 proteins are the most diverse in length and amino acid identity among hydrozoans (Fig. 2c, d). In D. pumila, gastrulation is apolar and mainly proceeds via epithelization of the outer cells. This mode of gastrulation causes deformations of the embryonic surface and results in multiple concavities and indentations. Thus, at the midgastrula stage, multiple epithelized toroidal surfaces compose an embryonic surface. These deformations are smoothed out towards the end of gastrulation, when only several indentations are still visible. The last indentation tends to be located in the oral domain of the embryo. However, this last indentation is not homologous to a blastopore 39 . At the end of gastrulation, in situ hybridization revealed expression of DpBra1 in a unitary broad domain (Fig. 3a) which did not overlap with any specific region within the gastrula stage embryo. Signal was visualized both in the ectoderm and the endoderm (Fig. 3b, c). At the preplanula stage, expression signal was detected in discrete patches both in the ectoderm and the endoderm mostly at the oral end of the embryo (Fig. 3d, e). In the early planula, we observed DpBra1 expression in the oral third of the larva (Fig. 3f, g). In the mature planula, DpBra1 was expressed in the oral half of the larva (Fig. 3h, i). In the ectoderm, we observed two domains of DpBra1 expression. In the tip, biased towards the oral pole expression signal was visualized in apical domains of ectodermal cells. Also, DpBra1 expression was visible in scattered ectodermal cells in the middle of the larva. In the endoderm, expression was present in only a few cells at the oral end (Fig. 3j). Figure 3k represents expression patterns of DpBra1 during development.  www.nature.com/scientificreports/ DpBra2 expression was detected both in the ectoderm and the endoderm at the end of gastrulation and also forms a broad domain (Fig. 4a, b). In the preplanula/early planula, DpBra2 expression was observed at the oral half of the embryo with more prominent signal at the oral pole (Fig. 4c). Transcripts concentrated in the perinuclear cytoplasm of ectodermal cells (Fig. 4d). In the mature planula, we observed DpBra2 expression in oral third of the larva with a bias towards the pole (Fig. 4e, f). DpBra2 RNA was visualized in apical domains of ectodermal cells (Fig. 4g). Expression was also present in single endodermal cells (Fig. 4g) as in case of DpBra1 (Fig. 3j). Figure 4h represents expression patterns of DpBra2 during development.
DpBra3 was expressed in a broad domain at the end of gastrulation (Fig. 5a) in a pattern similar to those of DpBra1 and DpBra2. However, the expression pattern of DpBra3 differed drastically at later developmental stages. At the preplanula stage, DpBra3 signal was visible as a belt in the center of the oral-aboral axis (Fig. 5b,  c). As the development proceeded, the expression area expanded to cover the central part of the early planula (Fig. 5d). Longitudinal sections (Fig. 5e, f) revealed that transcripts are present mostly in basal domains of scattered ectodermal cells. As the planula elongated, expression continued in the middle part of the larva in discrete ectodermal cells (Fig. 5g, h). Weak signal also appeared in the aboral endoderm (Fig. 5g, h). Longitudinal sections of the mature planula clearly demonstrated that bottle-like (Fig. 5i) and triangular (Fig. 5j) bodies of DpBra3expressing cells were located directly above the basal lamina or are between endoderm and ectoderm. The latter probably migrate towards the ectoderm from endoderm (Fig. 5h). Figure 5k represents expression patterns of DpBra3 during development.
Further, we examined expression patterns of the three Brachyury genes in the colony shoots of D. pumila. D. pumila forms monopodially growing colonies possessing biradial symmetry. Shoots of the colony are composed of repetitive modules. Each module consists of a fragment of the shoot in the center and two hydrants on the sides (Fig. 6a). New modules are formed on the top due to the repeating morphogenetic cycle in the specific organ-the shoot growth tip (Fig. 6b). Stage 1 represents the state when the morphogenetic cycle has not started yet (Fig. 6b1). At stage 2, the growth starts with the apical surface of the tip curving up (Fig. 6b2). At stage 3, the  www.nature.com/scientificreports/ tip elongates and takes a hemispherical shape (Fig. 6b3). At stage 4, the growth tip is dividing into the central and two lateral parts (Fig. 6b4). Lateral primordia further differentiate into hydrants, while the central part will become the new shoot growth tip (Fig. 6b1*). We analyzed expression patterns of three Brachyury genes in shoot growing tips on stages 2 and 4 of the morphogenetic cycle and in fully formed differentiated hydrants. DpBra1 and DpBra2 expression was detected in the apical ectoderm of the growth tip (Fig. 6c). DpBra1 is expressed in the central part of the apex at stage 2 and uniformly at stage 4. DpBra2 expression was observed in two domains at opposite sides of the apex at stages 2 and 4. Thus, DpBra1 and DpBra2 expression domains do not overlap at stage 2, but are co-expressed at stage 4. DpBra3 expression was not detected in the shoot growth tip.  www.nature.com/scientificreports/ Expression of three Brachyury genes was observed in the hypostome of the hydranth (Fig. 6d). A longitudinal section through the center of the hydranth revealed that DpBra1 was expressed both in the ecto-and the endoderm. DpBra2 signal was clearly visible in the endoderm of the hypostome, while the presence of a signal in the ectoderm is unclear. Unfortunately, DpBra3 signal was too weak for the fine examination, but seems to be expressed in the ectoderm (when viewed from the surface).
Brachyury genes are differently regulated by the cWnt signaling in D. pumila. It was shown previously in Hydra and C. hemisphaerica that two hydrozoan Brachyury genes, Brachyury1 and Brachyury2, are regulated by cWnt signaling 13,33,40,41 . However, it is unknown if Brachyury3 is still a cWnt-dependent gene after the duplication event. We assayed the dependence of three Brachyury genes on the cWnt pathway in D. pumila. We treated embryos at the gastrula stage with different concentrations of pharmacological agents to modulate the cWnt pathway, cultivated them until planula stage, and examined then expression patterns of three Brachyury genes in planula larvae of D. pumila (Fig. 7). Azakenpaullone (Azk) activates cWnt signaling and iCRT14 inhibits it [42][43][44] . It was shown in a previous study that hyper-activation of cWnt signaling results in the enlargement of larval oral domain, while its inhibition leads to reduction of oral domain in D. pumila 39 .
DMSO-treated (control) larvae had normal morphology and expression patterns of three Brachyury genes (Fig. 7). Treatments with the increasing concentrations of Azk resulted in the gradual expansion of DpBra1 and DpBra2 expression domains. After 2.5 μM Azk treatment, DpBra1 and DpBra2 expression signals were observed in the entire larva except the aboral-most region. The number of endodermal DpBra1-and Bra2-positive cells also increased. Vice versa, gradual inhibition of the cWnt signaling with iCRT14 led to the decrease of DpBra1 and DpBra2 expression domains in area (Fig. 7).
Strikingly, overactivation of the cWnt signaling did not lead to the expansion of DpBra3 expression domain. The belt of DpBra3-expressing cells shifted in more and more aboral positions, vacating the central domain. In the result of 2.5 μM Azk treatment, several DpBra3-expressing cells were detected at the aboral pole of the larva (Fig. 7: arrowhead). Inhibition of the cWnt signaling did not notably change the DpBra3 expression domain (Fig. 7).

D. pumila Brachyury genes differently regulate tissue differentiation in the animal cap assay.
To uncover functional differences of three D. pumila Brachyury genes, we employed the Xenopus laevis animal cap assay system. Using this assay, we surveyed DpBra1, DpBra2, and DpBra3 for their ability to affect cell fates of naive Xenopus animal cap cells. It is known, that untreated animal caps differentiate into epidermal tissue 45 , but the injection with Xenopus Bra or Hydra Bra1 mRNA promotes mesoderm specification, and Hydra Bra2 mRNA shows neural-inducing activity 33,46 . We injected capped mRNAs encoding DpBra1, DpBra2, or DpBra3 into the animal region of two-to four-cell stage embryos (~ 1 ng per embryo), dissected the animal caps www.nature.com/scientificreports/ at the blastula stage (stage 8), cultured them until control embryos reached late neurula stage (stage 18), and examined marker gene expression in these caps using conventional or quantitative RT-PCR (qRT-PCR). Uninjected animal caps were used as a control group. DpBra1 significantly (P < 0.0001) induced the expression of mesodermal marker gene actc1.L (muscle actin) 47 (Fig. 8a) as well the expression of another mesodermal marker gene, myod.S 47 . Since myod.S expression was too low for reliable quantification in the control group using qRT-PCR, we used gel electrophoresis to show the induction (Fig. 8c). DpBra1 did not affect the expression of neural marker gene tubb2b.S 47 (Fig. 8b) while DpBra2 and DpBra3 did not affect the expression of neuronal and mesodermal marker genes (Fig. 8).

Discussion
Gene duplications facilitate the evolution of regulatory genes, driving expansion of families of signaling molecules and transcription factors 48,49 . Duplication event produces two copies (paralogs) of a gene both of which are orthologous to the "parent" gene. Subsequent evolution may lead to divergence: one copy retains high similarity to its orthologs while the second undergoes structural changes and may be designated as a daughter or child copy 50 . In the present study, three paralogous genes of Brachyury transcription factor were found in hydrozoan lineage. The phylogenetic reconstruction suggests the first duplication event occurred before the hydrozoan clade branched off (in the common ancestor of meduzosoans or earlier). The "daughter" copy of the first duplication (Brachyury2/3) experienced then an additional duplication event in the common ancestor of hydrozoans (Fig. 1).
There are three main scenarios after gene duplication. First outcome leads to loss of function of duplicated copy which becomes a pseudogene or is lost. Two other outcomes are subfunctionalization and neofunctionalization 51,52 . In case of subfunctionalization, two duplicated copies share the original function of the ancestral gene, and both are required to preserve the entire ancestral function 51,53 . Overlapping expression domains of two duplicates often reflect occurred subfunctionalization 51 . In case of neofunctionalization, one copy retains its ancestral functions, and the other one is free to gain a novel function, since it is relaxed from the selective pressure 53 , and often acquires the expression domain different from the ancestral gene 54,55 . The acquisition of novel functions by regulatory genes plays a key role in diversification of developmental pathways and body plans in metazoans 56,57 . Importantly, the same duplicate can also display both features of sub-and neofunctionalization with regard to different functions 58 .
Evolution of Brachyury reveals signs of both sub-and neofunctionalization. In teleosts, expression of two Brachyury genes reveals the common chordate pattern 28,59 . Only simultaneous loss of both Brachyury genes recapitulates mouse homozygous brachyury mutant phenotype 4,28 , indicating subfunctionalization followed duplication of Brachyury in teleosts. On the contrary, neofunctionalization seems to follow Brachyury duplication in X. laevis as XBra3 (tbxt2.L/tbxt2.S) is distinct from XBra and XBra2 (tbxt,L/tbxt,S) in function and spatio-temporal pattern of expression 60 . Within hydrozoans, outcomes of Brachyury duplication previously were studied in Hydra 33 . Though both Brachyury paralogs are expressed in the hypostome of Hydra, they evolved distinct coding sequences and diverged their functions. Authors posit, that Brachyury paralogs show a mixture www.nature.com/scientificreports/ of sub-and neofunctionalisation in Hydra 33 . However, it is unknown, if similar outcomes have shaped the roles of Brachyury genes in other cnidarians. Brachyury genes are involved in axial patterning also in other studied cnidarians 14 . In cnidarian species with polar gastrulation, where it is linked to the axial patterning and occurs in the oral region of the embryo 61 , Brachyury expression accompanies gastrulation morphogenetic movements 12,19,20,36 . In gastrulating C. hemisphaerica, Brachyury paralogs (ChBra1 and ChBra2) display overlapping patterns 19 and are important for progression of gastrulation 13 . In D. pumila during gastrulation, all three Brachyury paralogs are expressed in a broad domain (Figs. 3a, 4a, 5a). In contrast to the cnidarian species with polar gastrulation 12 , it is unlikely, that Brachyury genes provide demarcation of ecto-endoderm boundary in D. pumila, since germ layers specification is not associated with axial polarity and oral region in particular during gastrulation in this species 39 .
Predominantly oral expression of Brachyury genes continues throughout a cnidarian life cycle. Brachyury is expressed in pharynxes of anthozoan larvae 12,20 and in the oral ectoderm of hydrozoan larvae in C. hemisphaerica and D. pumila, where expression of Brachyury1 and Brachyury2 is detected 19 (Figs. 3h-j, 4e-g). Surprisingly, Brachyury3 orthologs are not associated with oral tissues in hydrozoan larvae. Previously examined Brachyury ortholog of P. carnea which is expressed in the aboral ectoderm 36 clusters according to our analysis with Brachyury3 group (Fig. 1). In D. pumila, Brachyury3 displays expression in discrete triangular and bottle-like ectodermal cells (Fig. 5g-j), morphologically similar to sensory cells of cnidarian nervous system 62 , though Brachyury is known to act as neural repressor in the anthozoan Nematostella 14 . Moreover, DpBra3 seems not to be a cWnt-dependent which is in contrast to reported Wnt-dependence of Brachyury expression in cnidarians and bilaterians 14,63 including DpBra1 and DpBra2 (Fig. 6). Brachury3 orthologs display strong diversity in length and amino acid identity among hydrozoan Brachyury gene family (Fig. 2c, d). Differences in protein sequences, regulation, and expression domains suggest that the newly derived Brachyury3 diverged to different functions in hydrozoan species.
In the hydranth of D. pumila, DpBra3 expression pattern does not drastically differ from patterns of DpBra1 and DpBra2 and is in line with previous studies 33,36,64,65 . All three Brachyury paralogs were detected in hypostome of hydranth with overlapping patterns (Fig. 6d). Likely, the original function of Brachyury in a hydrozoan hydranth is associated with a specification of an oral domain (a hypostome) as a whole. As in Hydra 33 , overlapping expression domains of Brachyury paralogues in hypostomes of D. pumila hydranths suggest occurred subfunctionalization 51 .
In the shoot growing tip of D. pumila colonies 66 , DpBra1 and DpBra2 are strongly expressed in its apical ectoderm (Fig. 6c), which could be considered a derivative of the larval oral domain. Association of DpBra2 expression with the formation of hydranth primordia indicates its novel function in the hydranth primordia in D. pumila.
The specificity of Brachyury function is mostly defined by the N-and C-terminal domains, but not by the central T-box 9,67 . In line with previous studies 33 and our data (Fig. 8), high functional conservation of hydrozoan Brachyury1 orthologs is consistent with high conservation of protein sequence (Fig. 2c). However, our data indicates the functional divergence of Brachyury2 and Brachyury3. In Xenopus animal cap assays, DpBra2 and DpBra3 did not cause increased expression of the mesodermal markers actc1.L or myod.S (Fig. 8a, c). In Brachyury2 and Brachyury3, N-and C-terminal domains show lesser amino acid identity to the ancestral gene and have lost ancestral C-terminal repression domain R1 (Fig. 2a-c). These differences in terminal domains could be responsible for the neo-and subfunctionalization of Brachyury2 and Brachyury3 in hydrozoans, even though it was suggested that they occur mainly due to mutations in regulatory sequences, rather than mutations in the coding sequence 68 .
Taken together, our data indicate two duplication events of Brachyury in cnidarians. Brachyury1 is the most conservative duplicate, both on the functional and sequence levels. In studied hydrozoans and in D. pumila in particular, it is supposed to preserve its ancestral function as a crucial component of axis formation and patterning. Hydrozoan Brachyury 2 and Brachyury 3 reveal features of sub-and neofunctionalization. Brachyury3 however, displays strong divergence in sequence and functions among hydrozoans. Our data on Brachyury support the model of an indistinct border between sub-and neofunctionalization and complex outcomes for duplicated genes 58 , and provides a promising model for studies on post-duplication scenarios. Control embryos were exposed to 0.1% DMSO in FSW. Working solutions were refreshed daily. Incubation was performed in the dark.

Methods
Data sources and transcriptome assembly. To analyse phylogenetic relationships within the brachyury gene family, we surveyed 28 metazoan species. Gene sequences were obtained from several sources (Supplementary Information Table S1). Bilaterian, ctenophore, placozoan and anthozoan sequences were obtained  Information Table S1). We also used sequences of bilaterian, ctenophore, placozoan and anthozoan Brachyury genes. Nucleotide sequences with no corresponding protein sequence in the NCBI database were translated using Transdecoder v5.5.0. The search of T-boxes in analyzed sequences was performed with NCBI Conserved Domain Search tool. Amino acid sequence alignments and phylogenetic analysis were performed with MUSCLE algorithm in MUSCLE software (v3.8.31) 76 . Sequences of Tbx genes were selected as an outgroup. Sequence alignments were trimmed by removing poorly aligned regions using TrimAL tool, v.1.2rev59 77 . A heuristic approach "automated1" was used to select the best automatic method to trim our alignments. Trimming was performed without manual adjustment. Phylogenetic analysis was performed with Maximum Likelihood using IQTree v.2.0-rc2 software 78 . The JTT + R5 model was found to be optimal. To assess branch supports, bootstrap values were calculated running 1000 replicates using ultrafast bootstrap (UFBoot) 79 . Trees were visualized in FigTree v1.4.4 software. Obtained phylogenetic trees were processed with Adobe Illustrator CC. No corrections were made to the tree topology and the branch lengths.
We searched for Brachyury transcripts in gene models of Hydra 2.0 and HydraAEP genome assemblies using phylogenetically-informed annotation pipeline PIA3 80 . PIA3 pipeline is modified from PIA2 81 and is available on GitHub 82 . Modifications allowed us to automatically retrieve T-box protein class information.
To analyze functional domains of the hydrozoan Brachyury, selected protein sequences were scanned against Pfam hidden Markov model (HMM) database using hmmscan of HmmerWeb v.2.41.1 83 . Identification of the conserved R1 domain within the hydrozoan Brachyury was carried out using ClustalW sequence alignment service 84 with the R1 domain in HyBra1 33 as a query. The domain architecture of proteins was visualized using Pfam 85 . Multiple sequence alignment and calculation of the identity matrix of hydrozoan Brachyury proteins and T-boxes of D. pumila Brachyury were conducted using ClustalW with default settings and shaded using BOXSHADE 3.21.
D. pumila genes isolation, PCR, and antisense RNA probe synthesis. cDNA expression library was prepared by the SMART approach from total embryonic RNA with a Mint cDNA synthesis kit (Evrogen, Russia). cDNA gene fragments were isolated from the library by PCR with gene-specific primers (see Table 1). Primers were designed based upon sequences obtained from the sequenced transcriptome (Illumina) of D. pumila 35 . Amplified fragments were cloned into the pAL-TA vector (Evrogen, Russia). Digoxygenine-labeled antisense RNA probes were generated from gene fragments, which were amplified from plasmids with D. pumila genes.
In situ hybridization. The in situ hybridization protocol was performed as previously described in 66 for D. pumila shoots and hydranths and in 39 for D. pumila embryos. An urea-based in situ hybridization method was used for the hydranths 86 .
Shoots were fixed with 0.2% glutaraldehyde/4%formaldehyde in FSW for 1 min and then for an additional hour with 4% formaldehyde in FSW. Samples were washed with PTw (1× PBS with 0.1% Tween 20) thrice and stored in 100% methanol no more than overnight at −20 °C until hybridization. Embryos were fixed with 4% paraformaldehyde in FSW overnight at + 4 °C, rinsed with PBS, and stored at −20 °C in 100% methanol until hybridization. Table 1. PCR and qPCR primers used in this study.